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Abstract 

The aim of the article is to investigate the relative dispersion proper- 
ties of the Well Mixed class of Lagrangian Stochastic Models. Dimensional 
analysis shows that given a model in the class, its properties depend solely 
on a non-dimensional parameter, which measures the relative weight of 
Lagrangian-to-Eulcrian scales. This parameter is formulated in terms of 
Kolmogorov constants, and model properties are then studied by modify- 
ing its value in a range that contains the experimental variability. Large 
variations are found for the quantity g* = 2gCg 1 , where g is the Richard- 
son constant, and for the duration of the t 3 regime. Asymptotic analysis 
of model behaviour clarifies some inconsistencies in the literature and ex- 
cludes the Ornstein-Uhlenbeck process from being considered a reliable 
model for relative dispersion. 



1 Introduction 

Relative dispersion is a process that depends on the combination of the Eulerian 
and Lagrangian properties of turbulence. If particle separation falls in the iner- 
tial subrange, the Eulerian spatial structu re affects the dispersion, w hich can be 
regarded as a Lagrangian these property l|Monin and Yagloml |l975). The com- 
bination of properti es requires that both descriptions be considered (see e.g., 
Ifioffetta etallllfl^. 

Lagrangian Stochastic Modelling (LSM) is one turbulence representation 
that naturally combines Eulerian spatial structure and Lagrangian temporal cor- 
relation. In fact, as formulated bv lThomsonl l)l990l) using the Well Mixed Condi- 
tion (WMC), Lagrangian and Eulerian statistics are accounted for through the 
second order Lagrangian structure function and the probability density function 
(pdf ) of Eulerian velocity. Several studies prove that this approach leads to the 
qualitative repr oduction of the main properties, a s expected from the Richard- 
son theory (see iThomsonl ll99(A iRevnoldsl Il999t ISawfordl 120011 among oth- 
ers). Furthermore, recent experimental studies seem to confirm the validity of 
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the Markovianity ass umption for the velocity JPorta et all Effllt iRenner eTall 
l200lUMordant et allEonl . 

However, t he intrinsic non-uniqueness of the WMC formulation (see, e.g., 
[ SawfordL Il999l) and the indeterminacy of the Kolmogorov constants (see, e.g., 
brccnivasan, FflQaiAnfossi et a,lll200(l for reviews) do not allow for a completely 
reliable representation of the process. In particular, the value of the Richardson 
consta nt predicted by previous studies is not uniquely determined (see, among 
others, l^homsorlll990l:lKurbanmuradovlll997tlBorgas and SawfordLll994HRevnoldsl 



1999|. Whether this indetermination is a result of the different formulation of 
models, or of the different values of the parameters adopted is still unclear, and 
no systematic studies have been performed so far. 

It is worth noting that, even focusing attention only on the de pendence on 

the mo del constants produces significant variability. As an example. lBorgas and Sawfordl 
(1994,) present the variation of the Richardson constant value with the La- 
grangian Kolmogorov constant Co . 

The aim of this article is to investigate the general properties of models 
based on the WMC with regard to inertial subrange relative dispersion features. 
In Section [21 the properties of the WMC are evidenced through a dimensional 
analysis, while the limit for vanishing spatial correlation is studied in Section[3] 
Subsequently a model formulation is discussed in Section^ and results analysed 
in Sectional 



2 The non-dimensional form of the well mixed 
condition 

Following the logical development of iThomsonl l|l987l) . Frhomsonl l|l99(ll) (here- 
inafter T90) extended the method for the selection of single particle Lagrangian 
Stochastic Models to models for the evolution of particle pair statistics. In the 
latter models, the state of a particle pair is represented by the joint vector of 
position and velocity (x, u) = (x' 1 ', i6 2 \ u' 1 ', u^ 2 )), where the upper index de- 
notes the particle, whose evolution is given by the set of Langevin type equations 
(LE) (with implied summation over repeated indices): 

f dxi = Ui dt . . 

\ dui = a t (x, u, t ) dt + bij (x, t) dWj (t) , [ ' 

where i,j — 1..6. The coefficients a an d b are de t ermin ed, as usual, through 
the well known Well Mixed Condition l)ThomsonL Il987|) and the consistency 
with the inertial subrange scaling, respectively. Further details are not given 
here, in that they ar e well established and widely used in the literature (see, 
e.g., ISawfordl l20QlL for a review). The only remark we would make is that, 
although lThomsonl l|1987l) himself studied this alternative, the tensor bij cannot 
be dependent on u in order to allow Eq. 0} to describe a physical ly meaningful 
process. In fact, as shown, for instance, bv Ivan Kampenl l|l98lj) . the Ito and 
Stratonovich calculus give different results when bij — 6,y(u). In particular, the 
WMC would not have a unique definition. Thus, from now on, b. h j — \j C^sSij , 
i,j = 1.. 6 will be u s ed ac cording to the usual scaling of Lagrangian structure 
function llThornsonL lT987l). where e is the mean dissipation rate of turbulent 
kinetic energy. 
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It should be remembered here that the WMC is satisfi ed by constra ining; 
the Fokker-Planck equation associated to Eq. (see, e.g., iGardinerl Il990fl to 
be consistent with the Eulerian probability density function of the flow. In the 
case of particle pairs the considered pdf is the one-time, two-point joint pdf of 
and u^, i = 1, 2, accounting for the spatial structure of the turbulent flow 
considered. The open question about the non-u niqueness of the solution in more 
than one dimension (see, e.g., ISawfordll99<ih is not addressed here. However, 
the following analysis will show that the problem studied is independent of the 
particular solution selected. 

In order to highlight the effect of turbulence features on the model formula- 
tion, characteristic scales for particle pair motion must be identified. Because 
the process of relativ e dispersion has to deal wit h both Eulerian and Lagrangian 
properties (see, e.g., iMonin and YaglomLll975l p. 540), such scales can be de- 
fined by considering the second order Eulerian and Lagrangian structure func- 
tions, i.e., 

(Av 2 ) ~ C K (eAr) 2 / 3 (2) 

for Eulerian velocit y v for a separation Ar = ||Ar||, according to the standard 
IKolmogorovl ljl94ll) theory (hereinafter K41), and 

(Au 2 ) ~ C (et) (3) 

for Lagrangian velocity u (see, e.g., IMonin and Yagloml Il975)) . where Av = 
||v(r + Ar) — v(r)|| and Au = ||u(t + At) — u(i)||. A length scale A can be 
defined in the Eulerian frame, so that in the inertial subrange (namely, for 
rj <C T <C A where r\ is the Kolmogorov microscale) the structure function for 
each component may be written as 

(Av 2 ) = 2a 2 (^ (4) 



where a = -\/||v|| 2 /3, which together with Eq. J2J) provides a definition for A. 

A Lagrangian time scale r can be defined in a similar way using Eq. J3J and 
the Lagrangian version of Eq. (J3J. Thus, for r, « i < t, one has 



(Au 2 )=2a 2 - (5) 

T 

from which one can retrieve the known relationship 

suggested bv lTennekesI l|l982() . It should be observed that scales for the inertial 
subrange, at variance with their integral version, can be defined independently 
of non- homogeneity or unsteadiness, provided that the scales of such variations 
are sufficiently large to allow an inertial subrange to be identified. As far as the 
velocity is concerned, a can be recognised as the appropriate scale of turbulent 
fluctuations in both descriptions. 

The quantities er, A and r can then be used respectively to make velocity 
Ui, position Xi and time t non-dimensional. They also form a non-dimensional 
parameter 

or C«_ 

' A V2Q,' [> 
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the last equality being based on the combination of Eqs. and |JSJ with Q 
and JSJ). The parameter (3 can be recognised as a version of the well known 
Lagrangian-to-Eulerian scale ratio. The approach adopted here evidences its 
connection to fundamental constants of the K41 theory. 
In non-dimensional form, Eq. JQ) reads 



dxi = (3ui dt 

aui = a, dt + y/2dWi{t). 



(8) 



where, with a change of notation with respect to Eq. all the quantities 
involved are without physical dimensions. 
The associated Fokker-Planck equation is 

dpL | p u dp L | da^L = d 2 p L 
dt dxi dui duidui 

where pl is the pdf of the Lagrangian process described by Eq. (JHJ for some 
initial conditions. Using the WMC, a can be written as 



d\np E 

dui 

where 



(10) 



iPE dp E a dp E 

= —*r-P u i-azr ( n ) 



dui dt dxi 

and pe is the Eulerian one-time, two-point joint pdf of x and u. 

An advantage of this choice of scales emerges clearly in Eq. JHJ. It shows that, 
given a Eulerian pdf, once the non-uniqueness problem is solved by selecting 
a suitabl e solution to E q. (|1U[) . or applying a further physical constraint to 
Eq. (|llfl ( Sawfordl ll999). any solution of Eq. JOJ will depend on one parameter 
only, namely on the Lagrangian-to-Eulerian scale ratio. It can also be observed 
that this dependence is completely accounted for by the non-homogeneity term, 
which is an intrinsic property of the particle pair dispersion process in spatially 
structured velocity fields. 

In looking for the universal properties of pair-dispersion in the inertial sub- 
range, it is useful to rewrite the Richardson t 3 law in non-dimensional form, 
i.e., Ax 2 = g*f3 2 t 3 where g* = 2g/Co. In this form, the numerical value of the 
"normalised" Richardson constant g* depends on (3 only. This dependence is 
investigated in the following Sections to highlight the intrinsic properties of the 
LSM. 



3 The spatial decorrelation limit 

In the limit (3 — > oo, corresponding to a vanishing Eulerian correlation scale, the 
non-dimensionalisation defined in the previous section fails to apply. However, 
in this limit, the WMC solution can be proven to reduce to an homogeneous 
process (see Appendix). In particular, selecting a Gaussian pdf will give the 
Ornstein-Uhlenbeck (OU) process. It is worth noting that the OU process has 
sometime s been used to describe Lagrangian velocity in turbulent flows, for in- 
stance bv lGiffordl 1)1982) . who pioneered the stochastic approach to atmospheric 
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dispersion. The Novikov ||1963|) model and the NGLS model l|Thomsonl Fl990. 
p. 124) are simple applications of this concept. 

Adopting the choices made in the previous Section, but using the spatial 
scale defined by tct rather than the vanishing A as a length scale, the OU 
process equivalent to Eq. (jSJ is described by the non-dimensional set of linear 
LE 

dxi=Uidt , . 

dui = -Midi + a/2 dWi { ' 

where i — 1..6. The equations for the relative quantities (Aw^, Axi) can be 
obtained from the difference between quantities relative to the first (i = 1, 2, 3) 
and second (i = 4, 5, 6) particles. The resulting set of equations reads 

dAxi = Amdt . . 

dAu % = ~Au t dt + 2dW, [ ' 

where i = 1..3. 

Equation 113(1 can be solved analytically for correlation functions and vari- 
anc es (see e . g., \ Gardinerl [19901. Some basic results are summarised below (see 
also lOiffordLTl 9821). 

The second order moment of velocity difference turns out to be an exponen- 
tial function dependent on the time interval only 

((Am - Au m ) 2 ) = (Auli) cxp (-t) . (14) 

By integrating Eq. (|14fl . the displacement variance for a single component is 

((Axi - Ax 0% ) 2 ) = ({Au 2 0i ) - 2)(1 - exp (-t)) + At - 4(1 - exp (-*)) . (15) 

For short times (but expanding Eq. I|15|) to the third power of i), it turns out 
that 

((Aa* - Ax 0t ) 2 ) c (Au 2 0i )t 2 + (1 - (Au 2 0i )\ t 3 . (16) 

From Eq. i|16[l h can be observed that, when initial relative velocity Auoi is 
distributed in equilibrium with Eulerian turbule nce (i.e., ( Au 2 ^) = 2), a t 2 
regime takes place with a negative t 3 correction l|Hund . ll985|) . On the other 



hand, if (Au^) — the ballistic regi me displays a t 3 growth with a coeffic i ent 4, 
i.e., 2Cq for the dimension al version l|Novikovt ll963HMonin and Yagloml Il975t 
iBorgas and Sawfordlll99i]) . 



4 Model formulation and numerical simulations 

In order to proceed with the analysis of the dependence of model features on 
parameter (3, we select as a possible solution to Eq. I(l()|l . the expression given 
by T90 ( his eq. 1 8) for Gaussian pdf. The spatial structure is accounted for 
using the iDurbinl 0980) formula for longitudinal velocity correlation, which is 
compatible with the 2 /3 scaling law in the inertial subrange. Although this form 
is known not to satisfy completely the inertial subrange requirements (it pre- 
scribes a Gaussian distribution for Eulerian velocity differences, while inertial 
subrange r equires a non-zero skewnes s), it has been su ccessfully used i n ba- 
sic studies l|Borgas and Sawfordl Il994[) and applications llRevnoldsi ll9991. and 
provides a useful test case for studying the results shown above. 
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The stochastic model is formulated for the variable (x, u) rather than for the 
variable (Ax/^Au/v^) as in Thomson's original formulation. In the present 
case, assuming homogeneous and isotropic turbulence, the covariance matrix 
V(x) of the Eulerian pdf is expressed by 

V [ ^(x) 1 ) 
where X is the identity matrix and 

< P1 ' P2) (x) = K (pi) u? 2) ) (18) 

where pi,j»2 = 1, 2 (pi ^ p%) are the particles indices. The quantity {uf^u^^) = 
(ui(x.( pi * > )uj(x.( p2 " > )) is the two-point covariance matrix, which is e xpressed in 
terms of longitudinal and transverse functions F and G (see, e.g., iBatchelorl 
119531) as 

Kit = F(A x )Ax t Ax 3 + G{A x )5 i;j (19) 
where Ax = ||xW - x( 2 )||, 

and 

It goes without saying that ftg 1,Pa) = ft|f = ' P2) . As in lDurbinl jl98(fr . 
_F and G are computed from the parallel velocity correlation 

Ax 2 ^ 1/3 



which is K41 compliant for Ax <C 1. 

Using the above formulation, Eqs. © were solved numerically for a num- 
ber of trajectories large enough to provide reliable statistics for the relevant 
quantities. Particular attention was paid to the time-step-independence of the 
solution (details are not reported here). It was found that the time step strongly 
depends on [3 because large values of the parameter increase non-homogeneity, 
which requires greater acc uracy. Despite the widespread use of variable-time- 
step algorithms (see, e.g., lThomsonlll990l:ISchwere et all 12002^1 based, in par- 
ticular, on spatial derivatives, here a fixed time step short enough for time-step 
independence of the solution was used throughout the computation. 

Simulations were performed for two different initial conditions for velocity 
difference: i) the distributed case, where velocity differences are given according 
to the second-order Eulerian structure function and ii) the delta case ((Auf)o = 
0), where both particles of a pair arc released with the same velocity, which is 
normally distributed with variance 1. The two cases correspond to the limiting 
cases considered in Sectional The former describes "real" fluid particles, i.e., 
particles distributed like fluid at all times, while the latter represents, from the 
point of view of relative dispersion, marked particles leaving a "forced" source, 
where they were completely correlated (as for a jet). The initial condition for 
the spatial variable was Axq = 10~ 5 /3 for all simulations. It can be noted 
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that this corresponds to different positions in the inertial subrange for different 
simulations (Axq differs from case to case). However, the dimensional XAxo is 
chosen small enough to provide at least three decades of inertial subrange. 

The (3 parameter was varied in the range [10 -2 : 10 2 ], well beyond physically 
meaningful values . In fact, values reported in the literature range from OflO -1 ) 
to (X IO 1 ) llHinzel . Il95flt lHannaL Il98lt ISato and Yamamotol Il987t iKoeltzschl 
19990 with = 0(1) take n as a reference l|CorrsirTTl963(r This choice was 
made in order to infer asymptotic properties of the model. Note that, from 
a numerical point of view, different values of (3 were obtained by varying the 
length scale A, keeping a, r and Co fixed. In other words, with reference to 
Eq. JJJ, the variation of (3 was obtained by varying Ck- 



5 Results and discussion 

Figures ^ (from a to i) show the results of simulations for the two initial con- 
ditions and for different values of (3. The non-dimensional quantity (Ax 2 )f3~ 2 
is plotted against the non-dimensional time t. The OU analytical solutions 
{[3 = oo) are repor ted for referen ce. The general behavi our qualitatively filfills 
the expectations of lTavloil l|l92l|) and iRichardsonl ljl92(j|) . It presents an initial 
ballistic regime which differs for the two cases: the distributed case shows a 
t 2 , while the delta case presents a "false" 1 t 3 according to Eq. Hlfijl. After the 
ballistic regime there is a transition to an inertial range t 3 regime, which then 
becomes well established until a pure diffusive regime takes place. 

This generically correct behaviour merits further consideration. A "true" t 
is observed, which depends on the spatial flow structure and influences disper- 
sion properties. In particular, increasing [3 causes an increase in the normalised 
Richardson coefficient g* (Fig. |2J • It is worth noting that the "false" t 3 regime, 
according to the findings reported in Sect.|3| is not dependent on the structure , 
and therefore does not vary with [3. In fact, as pointed out bv lSawfordl i200lh . 
there should be a range where (Ax 2 )/3~ 2 = At 3 for t <C tn, to being the time 
at wh ich memory of the initial conditions is lost (see also lBorgas and Sawford 
Il99lh . It is clear now that this regime does not originate from any spatial struc- 
ture and is intrinsic to the solution with the delta initial condition, as explained 
by Eq. m- 

According to lMonin and Yagloml <|1975. p. 541), the "true" t 3 regime should 
be independent of the initial conditions. Thus, the starting point of this regime 
can be selected at the point where the solutions for the two cases coincide, as 
clearly occurs in Figs.^Ji to f. Therefore, the temporal extension of the t 3 regime 
is probably shorter than the one that could be estimated using intersections with 
the idealised ballistic regime, on the one hand, and with the diffusive regime, 
on the other. Note, however, that the extension of the inertial regime remains 
a decreasing function of (3, which asymptotically converges to zero. 

Another point of interest evident in Figs. a nd 2 is that the pres e nt re- 
sults do not agree with the theoretical findings of iBorgas and Sawfordl l|l99lj) 
(hereinafter BS91), although they compare well with the numerical results of 
IBorgas and Sawfordl l|l994h (hereinafter BS94). In fact, BS94 (their last figure) 
showed the results obtained by varying Co in their models. However, as shown 

1 In the sense that it is only a correction to the ballistic t 2 regime, which depends on the 
initial conditions and not on spatial structure. 
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Figure 1: Mean square separation normalised with f3 as a function of time 
for different values of (5. Thick lines represent results of present simulations, 
while thin lines are the analytical Ornstein-Uhlenbeck solutions (continuous: 
distributed case; dotted: delta case), a) /3 = 0.01, b) (i — 0.1, c) /3 = 0.2. 
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Figure 2: Normalized Richardson coefficient vs. f3. Pres ent results are repre- 
sented by ▲ connected with a line while V are taken from lBoreas and Sawfordl 
(1993). 



in Sect. (3 is the only parameter on which the model depends. Because of 
the constancy of e in BS94, the varation of Co corresponds to a variation of r, 
and hence (3. The results of BS94 for the implementation of the T90 model are 
reported in Fig. [2 and show a complete agreement with the present results. The 
values indicated do not satisfy the kinematic constraint g = 2Co — 7, where 7 
is a positive quantity, proposed by BS91, based on a double asymptotic expan- 
sion. It should be observed, however, that 7 is derived from kinematic features 
and depends on integrals of correlation functions. For vanishing correlation, one 
obtains 7 — > 0, suggesting that in BS91 the ballistic part of the OU process is 
an upper limit for dispersion in the Richardson regime. 

Nevertheless, this discrepancy can be explained as follows. From Eq. I|16|l 
it is clear that, at any time t < tdiff 5 where idiff is the time when pure dif- 
fusion takes place, the displacement variance for the OU process in the dis- 
tributed case is always larger than the displacement variance for the OU process 
in the delta case. The two cases represent the limit for any process based 
on the WMC for (3 — > 00, as shown in Sect. 03 In particular, focusing at- 
tention on the time range between the ballistic and diffusive regimes, it can 
be observed that lim^oo T90(distributed) = OU (distributed) and, further- 
more, T90(distributed) < OU (distributed) for any finite (3. It can be con- 
cluded that a (3' must exist for which T90(distributed) — OU (delta), and 
T90(distributed) > OU(delta) for (3 > j3' , in disagreement with BS91. Re- 
calling that the OU process is the limit for any WMC process with Gaussian 
Pe, this result can be considered to be applicable to more general kinematic 
properties, which should therefore depend on the ratio between Eulerian and 
Lagrangian scales. Thus, the limitation to g* in the BS91 derivation possibly 
derives from an implicit assumption concerning the spatial structure and/or the 
value of (3, which defines a range of applicability of the result. 

Proceeding further with the analysis, it can also be said that, because of the 
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existence of a time to after which the solution is not dependent on the initial 
conditions, it might be expected that T90(distributed) = T90(delta) for t > to- 
However, as the Lagrangian time increases with respect to (j _1 A, an increasing 
number of particle pairs reaches the end of the inertial range (Ax 3> 1) still 
remembering their initial conditions. This results in a range of (3 > 1 in which 
the delta solution never reaches the distributed solution before the onset of the 
diffusive regime. Therefore, it is not possible to define any g*. Nevertheless, for 
(3 > (3' there exists a range of t where T90(delta) > OU(deZta), which shows 
that T90(delta) converges to OV(delta) in a non-monotonic way. 

When, for f3 < 1, the expected independence on the initial conditions is 
recovered, it can be noted that to itself is a function of (3. Thus, the duration 
of the t 3 regime depends also (and mainly) on the starting time of the diffusive 
regime It is observed that decreasing (3 increases the time at which the diffusion 
regime becomes fully developed. 

6 Conclusions 

The dimensional analysis of the WMC, through the non-dimensionalisation of 
the Fokker- Planck equation has shown that only one parameter plays a role in 
the determination of two particle dispersion properties. This parameter is the 
Lagrangian-to-Eulerian scale ratio 0, which can be reliably defined in terms of 
inertial subrange constants. The dimensional analysis leads to the definition 
of a normalised Richardson constant g* whose scale is identified with Co, as 
suggested by the comparison of Lagrangian and Eulerian properties. Given a 
particular model, the numerical value of g* depends solely on the value of (3 
adopted. This also applies to the duration of the t 3 regime. 

Us ing the T90 formulation, it has been shown that the results of INovikovl 
(1963) are recovered for (3 oo, which means that in the model the spatial 
structure is negligible with respect to the Lagrangian time correlation. This 
limit corresponds to the OU process, whose general properties highlight that the 
observed t 3 growth is actually a correction to the ballistic regime t 2 . Moreover, 
because of the absence of any genuine t 3 regime, it is not possible to define any 
Richardson coefficient. This means that 2Co cannot be considered in general 
as the upper limit for g. Therefore there is no inconsistency in models that 
produce g > 2Co, as occurs in the present study and in BS94. 
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Appendix 

The stationary structure function of the second order, Eq. Q, can be generalized 
to an arbitrary integer order n, in non-dimensional terms as 

(Au n ) = (Au n ) c Ar hn , (A-l) 

where (-) c denotes Eulerian equilibrium statistics and, when n = 2, (Ait 2 ) c = 2. 
The inertial subrange and spatial decorrelation limit are recovered for h = 1/3 
and h = 0, respectively. 
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Considering the characteristic function pb(Aw; Ar) of the stationary Eule- 
rian pdf of velocity differences pe(Au; Ar) and using Eq. i|A-l(l . it turns out 
that 



p E (Aw; Ar) = ^(iAr' l Au;) n (Au n ) c (n!)- 1 = f(Ar h Aw), (A-2) 



n=0 



with i = a/— T- From Eq. l|A-2() it follows that 

MA-;Ar) = ^/(|^) , (A-3) 

where the factor Ar~ h conserves the normalization and, for the constant values 
h = 1/3,0, Eq. (|A-3j) defines the self similar regimes of the inertial subrange 
and the spatial decorrelation limit, respectively. 

Using the dimensional quantities Ar' = AAr and Au' — aAu for the particle 
separation and the velocity differences, respectively, for any finite Lagrangian 
correlation time r and particle separation Ar', the following identity holds 

lim Lp(Ar) = lim ip(Ar'/X) , (A-4) 

/3— >oc A^O 

where if is a generic continuous bounded function. Since continuity is required 
in the transition from the inertial subrange regime to the equilibrium, the scaling 
exponent h is assumed to be a monotonic decreasing function of Ar' /X. Thus 

lim A /l = 1 . (A-5) 

As observed in Sectional the only term affected by variations of (3 in Eq. (0 
is the non-homogeneous one. Therefore for any finite Ar' using Eq. I|A-3J) and 
Eq. HA-5(1 . it turns out that 



hm (3—- 

p^oc Or 



y (. h h ,( Au'a- 1 \ 



,2h+l 

which shows that the non-homogeneous term vanishes in this limit 
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